function agent_data_struct = natural_policy_gradient(agent_data_struct,...
                                 policy_derivatives)
    
    M = 1;
    % Check if converged
%     npg_struct.converged = has_converged_relatively(npg_struct.policy_gradient(:,end-1),...
%                                 npg_struct.policy_gradient(:,end), npg_struct.min_relative_change);
%     if npg_struct.converged
%         % converged already -> just return
%         return;
%     end
    observations = agent_data_struct.observations;
    actions = agent_data_struct.actions;
    rewards = agent_data_struct.rewards;
    theta = agent_data_struct.theta;
    % compute sufficient statistics
    if agent_data_struct.time_variant_baseline
        error('Not yet implemented');
    else
        [psi, F, g, phi, rbar] = sufficient_stats(policy_derivatives, agent_data_struct.gamma, observations, actions, rewards,theta);
    
        Q = (1+(phi'*((M*F - phi*phi')\theta)))/M;
        b=Q*(rbar - phi'*(F\g));
        gng = F\(g-phi*b);
        agent_data_struct.policy_gradient = [agent_data_struct.policy_gradient, gng];
%         if has_converged_relatively(agent_data_struct.gng, gng, agent_data_struct.min_relative_change)
%             agent_data_struct.policy_gradient = [agent_data_struct.policy_gradient, gng];
%             agent_data_struct.trials = 0;
%             agent_data_struct.converged = true;
%         end
    end
end

function [psi, F, g, phi, rbar] = sufficient_stats(policy_derivatives, gamma, observations, actions, rewards,theta)
%SUFFICIENT_STATS calculates the sufficient statistics for the natural
% log_policy_gradient_fun: function to calculate the gradient of the log
%                          of the policy 
% psi: policy derivatives
% F: Fisher matrix
% g: Vanilla gradient
% phi: Eligibility
% rbar: Average reward

    gamma_vect = [1, cumprod(repmat(gamma, 1, size(observations, 2)-1))];
    psi = policy_derivatives(theta, observations, actions);
    F = sum(psi, 2) * sum(psi, 2)';
    g = sum(psi,2) * (gamma_vect*rewards');
    phi = sum(psi, 2);
    rbar = (gamma_vect*rewards');
end
